/*

  Simpson's 1/3 rd and 3/8
Author:Vishwajeet Dusane
*/


#include <iostream.h>
#include <math.h>
#define TOLERANCE 0.0001
using namespace std;
float f1x(float x)
{
	return 1/(1 + (x*x));
};

float f2x(float r,float x)
{
	return sqrt(r*r - x*x);
};

class SimpsonMethod1/3Rule
{
	private : float fun_val[10000];
		  float x_val[10000];
		  float pi[100],h,lower_bound,upper_bound;
		  int n;
	public: SimpsonMethod1/3Rule();
                void integrate(int iterations,int choice);
		void generate_x_val();
		void generate_fun_val_by_f1x();  
		void generate_fun_val_by_f2x();  
		void show_pi_value(int iterations);
		int check_tolerance_value(int iterations);
		void setn();
		void set_bounds();
		void set_bounds(float radius);
	friend float f1x(float x);
};

class SimpsonMethod3/8Rule
{
	private : float fun_val[10000];
		  float x_val[10000];
		  float pi[100],h,lower_bound,upper_bound;
		  int n;
	public: SimpsonMethod3/8Rule();
                void integrate(int iterations,int choice);
		void generate_x_val();
		void generate_fun_val_by_f1x();  
		void generate_fun_val_by_f2x();  
		void show_pi_value(int iterations);
		int check_tolerance_value(int iterations);
		void setn();
		void set_bounds();
		void set_bounds(float radius);
	friend float f1x(float x);
};



SimpsonMethod1/3Rule::SimpsonMethod1/3Rule()
{
	n = 9;
}


void SimpsonMethod1/3Rule::generate_x_val()
{
	h = (float)(upper_bound - lower_bound)/n;
	x_val[0] = lower_bound;
	int i=1;
	for(i=1;i<n;i++)
		x_val[i] = lower_bound + (i*h);
}

void SimpsonMethod1/3Rule::generate_fun_val_by_f1x()
{
	int i;
	for(i=0;i<n;i++)
		fun_val[i] = f1x(x_val[i]);
}

void SimpsonMethod1/3Rule::generate_fun_val_by_f2x()
{
	int i;
	for(i=0;i<n;i++)
		fun_val[i] = f2x(upper_bound,x_val[i]);
}

void SimpsonMethod1/3Rule::integrate(int iterations,int choice)
{
	int i;
	pi[iterations] = (float)(h/3) * (fun_val[0] + fun_val[n-1]);
	for(i=1;i<n-1;i++)
	{
		if(i%2 == 0)
			pi[iterations]  += (float)(2*h/3)*fun_val[i];
		else
			pi[iterations]  += (float)(4*h/3)*fun_val[i];
	}		
	pi[iterations] *= 4;
	
	if(choice == 2)
		pi[iterations] /= (upper_bound * upper_bound);
}
void SimpsonMethod1/3Rule::show_pi_value(int iterations)
{
	cout << "N : " << n << "\t\tPI : " << pi[iterations] << endl;  		
}

int SimpsonMethod1/3Rule::check_tolerance_value(int iterations)
{
	float diff = pi[iterations] - pi[iterations - 1];
	if(diff < 0)
		diff *= -1;
	if(diff < TOLERANCE)		                                                                                                    return 0;
	else
		return 1;
}

void SimpsonMethod1/3Rule::setn()
{
	n *= 2;
}

void SimpsonMethod1/3Rule::set_bounds()
{
	lower_bound = 0;
	upper_bound = 1;
}

void SimpsonMethod1/3Rule::set_bounds(float radius)
{
	lower_bound = 0;
	upper_bound = radius;
}
/*******************************Functions for Simpsons 3/8th rule*************************************/
SimpsonMethod3/8Rule::SimpsonMethod3/8Rule()
{
	n = 8;
}


void SimpsonMethod3/8Rule::generate_x_val()
{
	h = (float)(upper_bound - lower_bound)/n;
	x_val[0] = lower_bound;
	int i=1;
	for(i=1;i<n;i++)
		x_val[i] = lower_bound + (i*h);
}

void SimpsonMethod3/8Rule::generate_fun_val_by_f1x()
{
	int i;
	for(i=0;i<n;i++)
		fun_val[i] = f1x(x_val[i]);
}

void SimpsonMethod3/8Rule::generate_fun_val_by_f2x()
{
	int i;
	for(i=0;i<n;i++)
		fun_val[i] = f2x(upper_bound,x_val[i]);
}
void SimpsonMethod3/8Rule::integrate(int iterations,int choice)
{
	int i;
	pi[iterations] = (float)(3*h/8) * (fun_val[0] + fun_val[n-1]);
	for(i=1;i<n-1;i++)
	{
		if(i%3 == 0)
			pi[iterations]  += (float)(3*h/8)*2*fun_val[i];
		else
			pi[iterations]  += (float)(3*h/8)*3*fun_val[i];
	}		
	pi[iterations] *= 4;
	if(choice == 2)
		pi[iterations] /= (upper_bound * upper_bound);
}
void SimpsonMethod3/8Rule::show_pi_value(int iterations)
{
	cout << "N : " << n << "\t\tPI : " << pi[iterations] << endl;  		
}

int SimpsonMethod3/8Rule::check_tolerance_value(int iterations)
{
	float diff = pi[iterations] - pi[iterations - 1];
	if(diff < 0)
		diff *= -1;
	if(diff < TOLERANCE)		                                                                                                    return 0;
	else
		return 1;
}

void SimpsonMethod3/8Rule::setn()
{
	n *= 2;
}

void SimpsonMethod3/8Rule::set_bounds()
{
	lower_bound = 0;
	upper_bound = 1;
}

void SimpsonMethod3/8Rule::set_bounds(float radius)
{
	lower_bound = 0;
	upper_bound = radius;
}
int main()
{
	Trapezoidal t1;
	SimpsonMethod1/3Rule s_13;
	SimpsonMethod3/8Rule s_38;
	int i=1,choice;
	float radius;
	cout << "****************************************************************************************" << endl;
	cout << "                            MENU" << endl;
	cout << "****************************************************************************************" << endl;
	cout << "1.Function is 1/(1+ x^2)" << endl;
	cout << "2.Function is sqrt(r^2 - x^2)" << endl;
	cout << "****************************************************************************************" << endl;
	cout << "Enter your choice : ";
	cin >> choice;
	if(choice == 1)
	{
		t1.set_bounds();
		s_13.set_bounds();
		s_38.set_bounds();
	}
	else if(choice == 2)
	{
		cout << "Enter the radius of circle : ";
		cin >> radius;
		t1.set_bounds(radius);
		s_13.set_bounds(radius);
		s_38.set_bounds(radius);
	}
	else
	{
		cout << "Bad Choice.Exiting the program." << endl;
		exit(0);
	}

	cout << "N indicating no. of intervals,value of PI generated using Trapezoidal Rule is : " << endl;
	while(1)
	{
		t1.generate_x_val();
		if(choice == 1)
			t1.generate_fun_val_by_f1x();
		else
			t1.generate_fun_val_by_f2x();
		
		t1.integrate(i,choice);
		
		t1.show_pi_value(i);
		
		if(i!=0)
		{
			if(!t1.check_tolerance_value(i))
				break;
		
		}
		t1.setn();
		i++;
	}

	cout << "N indicating no. of intervals,value of PI generated using Simpsons 1/3rd rule is : " << endl;
	while(1)
	{
		s_13.generate_x_val();
		if(choice == 1)
			s_13.generate_fun_val_by_f1x();
		else
			s_13.generate_fun_val_by_f2x();
		
		s_13.integrate(i,choice);
		
		s_13.show_pi_value(i);
		
		if(i!=0)
		{
			if(!s_13.check_tolerance_value(i))
				break;
		
		}
		s_13.setn();
		i++;
	}
	cout << "N indicating no. of intervals,value of PI generated using Simpsons 3/8th rule is : " << endl;
	while(1)
	{
		s_38.generate_x_val();
		if(choice == 1)
			s_38.generate_fun_val_by_f1x();
		else
			s_38.generate_fun_val_by_f2x();
		
		s_38.integrate(i,choice);
		
		s_38.show_pi_value(i);
		
		if(i!=0)
		{
			if(!s_38.check_tolerance_value(i))
				break;
		
		}
		s_38.setn();
		i++;
	}
	return 0;
}
